function varargout = View4D_BRIK(varargin)
% Viewer for 4D .BRIK file.
% - Allows for montages across both slice dimensions and time.
% - Also displays timeseries for selected voxel.
%
% Usage: 
%  View4D_BRIK(anat_img, brik_img, cellarray_of_other_funcdata, custom_legendnames)
%
% Required Inputs:
%  anat_img: Is your anatomical+tlrc.BRIK that serves as underlay.
%  brik_img: Is your the functional+tlrc.BRIK file that serves as overlay.
%  
% Optional Inputs (Can leave empty if not desired):
%  cellarray_of_other_funcdata:
%    - In the viewer, timeseries data of the selected voxel for the primary functional image
%      is plotted in a solid line. If you wish to overlay timeseries data from other functional
%      data onto the plot, you can specify the other datasets in the third argument.
%    - If left empty, only the primary functional image will be plotted.
%
%  custom_legendnames:
%    - Allows you to specify your own legend for the voxel timeseries plot.
%    - If left empty, dataset names will be used.
%
%
% Examples:
%  View4D_BRIK('anat+tlrc.BRIK', 'functional+tlrc.BRIK', [], [])
%
%  View4D_BRIK('anat+tlrc.BRIK', 'Cond1+tlrc.BRIK', ...
%      {'Cond2+tlrc.BRIK, ControlCond+tlrc.BRIK'}, {'Cond1', 'Cond2', Control'})
%
% Written by: Natasa Kovacevic
% Modified in areas marked with ** (Nov, 2013: M. Cheung)
% Last modified: Jan. 15, 2014


% Copyright (C) 2013-2014, Natasa Kovacevic
%
% This file is a part of the MEG & PLS Pipeline (MEGPLS). For more 
% details, see the documentation included with the software package.
%
% MEGPLS is free software: you can redistribute it and/or modify it under
% the terms of the GNU General Public License version 2 as published by 
% the Free Software Foundation. This program is distributed in the hope 
% that it will be useful, but WITHOUT ANY WARRANTY; without even the 
% implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. 
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License along
% with this program. If not, you can download the license here: 
% <http://www.gnu.org/licenses/old-licenses/gpl-2.0>.


% NK: note use set/get for fields that come prediscribed with objects (e.g. sliders) 
% and for fields that I added to handles. Use setappdata/getappdata for gui
% Begin initialization code - DO NOT EDIT
gui_Singleton = 0;
gui_State = struct('gui_Name',       mfilename, ...
   'gui_Singleton',  gui_Singleton, ...
   'gui_OpeningFcn', @View4D_BRIK_OpeningFcn, ...
   'gui_OutputFcn',  @View4D_BRIK_OutputFcn, ...
   'gui_LayoutFcn',  [] , ...
   'gui_Callback',   []);
if nargin && ischar(varargin{1})
   gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
   [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
   gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT


%--- Executes just before View4D_BRIK is made visible. ---%
%---------------------------------------------------------%
function View4D_BRIK_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to View4D_BRIK (see VARARGIN)

% Choose default command line output for View4D_BRIK
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% Check inputs:
if numel(varargin) < 4 | numel(varargin)>4
   disp('You must supply: anatomical underlay and functional 4D .BRIK file');
   disp('All files should be in AFNI format');
   disp('Usage: View4D_BSR(anat_img, brik_img, cellarray_of_otherdata, customlegendnames)');
   disp('Example: View4D_BSR(''anat+tlrc.BRIK'', ''functional+tlrc.BRIK'', [], [])');
   error('wrong number of inputs ');
end

% load anatomical BRIK file
Opt.format = 'vector';
[err,anat,Info,ErrMessage]=BrikLoad(varargin{1},Opt);
handles.anat = anat;

% load functional 4D .BRIK image:
Opt.format = 'vector';
[err,img,Info,ErrMessage]=BrikLoad(varargin{2},Opt);
BSR.img = img;
handles.dims = size(BSR.img);  % 4-dimensional (x, y, z, time)
if numel(size(BSR.img))==3, handles.dims(4) = 1; end
handles.BSR = BSR;
handles.Info = Info;

handles.MainFuncFile = varargin{2};
guidata(hObject, handles);

%** Add ability to overlay multiple datasets in voxel timeseries axes.
% load other .BRIK data to include in voxel timeseries:
if ~isempty(varargin{3})
    for l=1:length(varargin{3})
        Opt.format = 'vector';
        [err,OtherImgData{l},Info,ErrMessage]=BrikLoad(varargin{3}{l},Opt);
        
        % Check dims of other image data.
        if size(OtherImgData{l}) ~= size(BSR.img)
            error('Optional datasets must be same dimension as the main functional image.');
        end
    end
    handles.OtherImgData = OtherImgData;
    handles.OtherImgName = varargin{3};
else
    handles.OtherImgData = [];
    handles.OtherImgName = [];
end

%** Add ability to specify own legend in voxel-timeseries plot:
if ~isempty(varargin{4})
    handles.CustomLegend = varargin{4};
    if length(handles.CustomLegend) ~= length(handles.OtherImgData) + 1
        disp('Error: Custom legend does not contain correct number of inputs.');
        disp(' - Legend names should include both the main functional image and optional datasets.');
        error('Custom Legend has wrong number of inputs.');
    end
else
    handles.CustomLegend = [];
end
    

% Sets title to PLS LV info:
%tmp=strfind(Info.HISTORY_NOTE, '***');
% if ~isempty(tmp)
%     tmp=Info.HISTORY_NOTE(tmp:end);
%     loc=strfind(tmp,'\n');
%     finaltmp=tmp(1:(loc(1)-1));
%     set(handles.Title,'String',finaltmp);
% else

    %** Changed to just display full filepath as title
    set(handles.Title,'String',varargin{2});
%end

% set view and sliders
x = ceil(handles.dims(1)/2);
y = ceil(handles.dims(2)/2);
z = ceil(handles.dims(3)/2);

set(handles.AxialSlider,'Min',1,'Max',handles.dims(3),'SliderStep',[1/(handles.dims(3)-1) 1/(handles.dims(3)-1)],'Value',z);
set(handles.SagitalSlider,'Min',1,'Max',handles.dims(1),'SliderStep',[1/(handles.dims(1)-1) 1/(handles.dims(1)-1)],'Value',x);
set(handles.CoronalSlider,'Min',1,'Max',handles.dims(2),'SliderStep',[1/(handles.dims(2)-1) 1/(handles.dims(2)-1)],'Value',y);

setappdata(handles.Axialview,'slicenum',z);
setappdata(handles.Sagitalview,'slicenum',x);
setappdata(handles.Coronalview,'slicenum',y);

if numel(size(BSR.img))==3
   handles.dims(4) = 1;
   set(handles.TimeptSlider,'Min',1,'Max',1.001,'SliderStep',[1 1],'Value',1);
else
   set(handles.TimeptSlider,'Min',1,'Max',handles.dims(4),'SliderStep',[1/(handles.dims(4)-1) 1/(handles.dims(4)-1)],'Value',1);
end

% set up montage selection listbox and monatge viewing parameters
set(handles.SelectMontage,'String',{'None','Axial','Sagittal','Coronal'});
set(handles.MontageViewingTable,'data',[1 1 handles.dims(1); 1 1 handles.dims(2); 1 1 handles.dims(3); 1 1 3;]);
% start with first timept
set(handles.Timept,'String','1');

% set up bsr, tbsr, blend
UpdateBlend(handles);

% start using gui
UpdateSlices(x,y,z,handles);
UpdateTimept(1,handles);
guidata(hObject, handles); % Update handles structure


function varargout = View4D_BRIK_OutputFcn(hObject, eventdata, handles)
varargout{1} = handles.output;


function UpdateBlend(handles)
% set up bsr, tbsr, blend
% this function is normally called only in the begining and when user
% changes bsr thresholds -> because thhreshold is the only thing that
% alters the blend and the colormap
% clamp bsr image for better viewing

%** get clamp range from textbox:
%bsr_clamp = 5;
bsr_clamp = str2num(get(handles.TextboxClampRange, 'String'));
bsr = handles.BSR.img;
bsr(find(bsr>=bsr_clamp)) = bsr_clamp; bsr(find(bsr<=-bsr_clamp)) = -bsr_clamp;


% get bsr threshold
%** Added check to prevent threshold from exceeding max/min of data or clamp
posthresh = str2num(get(handles.BSR_posthresh,'String'));
negthresh = str2num(get(handles.BSR_negthresh,'String'));

maxvalue = getappdata(handles.BSR_posthresh,'d');
minvalue = getappdata(handles.BSR_posthresh,'c');

if posthresh > maxvalue | posthresh <= minvalue
    posthresh = maxvalue;
    set(handles.BSR_posthresh,'String', posthresh);
end
if negthresh < minvalue | negthresh >= maxvalue
    negthresh = minvalue;
    set(handles.BSR_negthresh,'String', negthresh);
end

% zero out bsr voxels below threshold
tbsr = bsr;
tbsr(find(tbsr<posthresh & tbsr>negthresh)) = 0;
orig_tbsr = tbsr;


% now map tbsr for colormap
d = max(max(max(max(bsr)))); c = min(min(min(min(bsr))));
if c==d, % that would happen only if tbsr is zero everywhere
   tbsr = ones(size(tbsr));
else
   tbsr = round(255*(tbsr-c)/(d-c)+1);
end

[cmap, ignore_pts] = set_colormap(d,c,posthresh, negthresh);
colormap(cmap);

% map grayscale image into midwindow in the colormap
img = handles.anat;
a = min(min(min(img)));
b = max(max(max(img)));
img = round(ignore_pts(1) + (ignore_pts(end)-ignore_pts(1))*(img-a)/(b-a));
img = min(max(img,ignore_pts(1)),ignore_pts(end));
img = repmat(img,[1 1 1 handles.dims(4)]);
blend = zeros([size(tbsr) 3]);

for dim=1:3
   utmp = reshape(cmap(img,dim),size(tbsr));
   otmp = reshape(cmap(tbsr,dim),size(tbsr));
   blend(:,:,:,:,dim)= (orig_tbsr==0).* utmp + (orig_tbsr~=0).* otmp;
end
% attach blend and colormap related infor to be used in other gui
% components
setappdata(handles.BSR_posthresh,'blend',blend);
setappdata(handles.BSR_posthresh,'cmap',cmap);
setappdata(handles.BSR_posthresh,'ignore_pts', ignore_pts);
setappdata(handles.BSR_posthresh,'c',c);
setappdata(handles.BSR_posthresh,'d',d);

function UpdateSlices(x,y,z,handles)
% get current bsr thresholds
posthresh = str2num(get(handles.BSR_posthresh,'String'));
negthresh = str2num(get(handles.BSR_negthresh,'String'));
bsr = handles.BSR.img;
% get  current time point and make sure that it makes sense
tpt = round(get(handles.TimeptSlider,'value'));

% get blend and colormap
blend_allt= getappdata(handles.BSR_posthresh,'blend');
blend = squeeze(blend_allt(:,:,:,tpt,:));
cmap = getappdata(handles.BSR_posthresh,'cmap');
colormap(cmap);

% AFNI IJK coordinates:
set(handles.Coord,'String',[num2str(x) '  ' num2str(y) '  ' num2str(z)]);

% AFNI indexing convention for AFNI viewer & functions starts at 0.
% Note: AFNI reads orientation in voxel-storage order, so LPI in AFNI is RAS (in spatial direction)
% Recall: RAS is orientation of coordinate systems of TLRC and MNI templates.
[err, XYZdic] = AFNI_Index2XYZcontinuous([x-1, y-1, z-1], handles.Info, 'LPI');
set(handles.MNI_coord, 'String', [num2str(XYZdic(1)) '  ' num2str(XYZdic(2)) '  ' num2str(XYZdic(3))]);

% udate slider positions
set(handles.AxialSlider,'Value',z);
set(handles.SagitalSlider,'Value',x);
set(handles.CoronalSlider,'Value',y);

setappdata(handles.Axialview,'slicenum',z);
setappdata(handles.Sagitalview,'slicenum',x);
setappdata(handles.Coronalview,'slicenum',y);

% update slice views
col = [0 0.8 0];
axes(handles.Axialview); % make this current axes
tmp=permute(squeeze(blend(:,:,z,:)),[2 1 3]);
image(tmp,'ButtonDownFcn',{@click_axial,handles});
set(handles.Axialview,'ydir','normal','xtick',[],'ytick',[]);
line([x x],[1 handles.dims(2)],'color',[0 1 0],'visible','on');
line([1 handles.dims(1)],[y y],'color',[0 1 0],'visible','on');
text(2,2,['z=' num2str(z)],'fontsize',12,'fontweight','bold','color',[0 1 0.3]);

axes(handles.Sagitalview); % make this current axes
image(permute(squeeze(blend(x,:,:,:)),[2 1 3]),'ButtonDownFcn',{@click_sagital,handles});
set(handles.Sagitalview,'ydir','normal','xtick',[],'ytick',[]);
line([y y],[1 handles.dims(3)],'color',[0 1 0],'visible','on');
line([1 handles.dims(2)],[z z],'color',[0 1 0],'visible','on');
text(2,2,['x=' num2str(x)],'color','g','fontsize',12,'fontweight','bold');
text(0,handles.dims(3)+2,num2str(bsr(x,y,z,tpt)),'fontsize',12,'fontweight','bold','color','c');

axes(handles.Coronalview); % make this current axes
image(permute(squeeze(blend(:,y,:,:)),[2 1 3]),'ButtonDownFcn',{@click_coronal,handles});
set(handles.Coronalview,'ydir','normal','xtick',[],'ytick',[]);
line([x x],[1 handles.dims(3)],'color',[0 1 0],'visible','on');
line([1 handles.dims(1)],[z z],'color',[0 1 0],'visible','on');
text(2,2,['y=' num2str(y)],'color','g','fontsize',12,'fontweight','bold');
text(0,handles.dims(3)+2,num2str(bsr(x,y,z,tpt)),'fontsize',12,'fontweight','bold','color','m');

% update BSR colorbar
axes(handles.colorbar_view);
imagesc([1:256]');  colormap(cmap);
ignore_pts = getappdata(handles.BSR_posthresh,'ignore_pts');

%**Commented out: Modified below to be more robust accomodating different datasets
%set(gca, 'ydir', 'normal', 'ytick',[1 ignore_pts(1) ignore_pts(end) 256], 'YTickLabel', {num2str(getappdata(handles.BSR_posthresh,'c')), num2str(negthresh), num2str(posthresh), num2str(getappdata(handles.BSR_posthresh,'d'))}, 'XTickLabel', []);
if ignore_pts(1) <= 1
    set(gca, 'ydir', 'normal', 'ytick',[1 ignore_pts(end) 256], 'YTickLabel', {num2str(getappdata(handles.BSR_posthresh,'c')), num2str(posthresh), num2str(getappdata(handles.BSR_posthresh,'d'))}, 'XTickLabel', []);
elseif ignore_pts(end) >= 256
    set(gca, 'ydir', 'normal', 'ytick',[1 ignore_pts(1) 256], 'YTickLabel', {num2str(getappdata(handles.BSR_posthresh,'c')), num2str(negthresh), num2str(getappdata(handles.BSR_posthresh,'d'))}, 'XTickLabel', []);
elseif ignore_pts(1) <= 1 && ignore_pts(end) >= 256
    set(gca, 'ydir', 'normal', 'ytick',[1 256], 'YTickLabel', {num2str(getappdata(handles.BSR_posthresh,'c')), num2str(getappdata(handles.BSR_posthresh,'d'))}, 'XTickLabel', []);
else
    set(gca, 'ydir', 'normal', 'ytick',[1 ignore_pts(1) ignore_pts(end) 256], 'YTickLabel', {num2str(getappdata(handles.BSR_posthresh,'c')), num2str(negthresh), num2str(posthresh), num2str(getappdata(handles.BSR_posthresh,'d'))}, 'XTickLabel', []);
end


% update BSR time series plot
axes(handles.BSR_timeseries);
Tseries = squeeze(handles.BSR.img(x,y,z,:));
if numel(Tseries) > 1
   at = min(Tseries); bt = max(Tseries);
   num_tpts = numel(Tseries);
   tpts = [1:num_tpts];
   plot(tpts,Tseries, 'LineStyle', '-', 'LineWidth', 2);
   hold on;
   if bt > at
      e = (bt-at)/20;
      set(gca,'ylim',[at-e bt+e],'xlim',[0 numel(Tseries)+1]);
   else
      set(gca,'ylim',[-1 1],'xlim',[0 numel(Tseries)+1]);
   end
   hold all;
end

%** Add ability to overlay multiple datasets in voxel timeseries axes.
% overlay other specified datasets onto timeseries plot
if ~isempty(handles.OtherImgData)
    num_OtherData = length(handles.OtherImgData);
    tmparray = zeros(num_OtherData, length(Tseries));
    for l=1:num_OtherData
        tmparray(l,:) = squeeze(handles.OtherImgData{l}(x,y,z,:));
    end
    tmparray = tmparray';
    plot(tmparray, 'LineStyle', '-', 'LineWidth', 1); 
    hold on;
    
    %** add ability for custom legendname
    if ~isempty(handles.CustomLegend)
        clickableLegend(handles.CustomLegend, 'interpreter', 'none');
    else
        [~, FuncFile, ~] = fileparts(handles.MainFuncFile);
        LegendData{1} = FuncFile;
        for f = 1:num_OtherData
            [~, OtherFile, ~] = fileparts(handles.OtherImgName{f});
            LegendData{f+1} = OtherFile;
        end
        clickableLegend(LegendData, 'interpreter', 'none');
    end
end

% horizontal dotted red lines indicate bsr thresholds
h=line([1 num_tpts], [posthresh posthresh]);
set(h, 'LineStyle', '--', 'Color', 'r');
h=line([1 num_tpts], [negthresh negthresh]);
set(h, 'LineStyle', '--', 'Color', 'r');
hold off;
title(['Value = ',num2str(Tseries(tpt))],'color','r','fontsize',14);
minval = min(min(Tseries),negthresh); maxval = max(max(Tseries),posthresh);
E = (maxval - minval)/10;
set(gca,'ylim',[minval-E maxval+E]);
% vertical gray line indicates current tpt
h = line([tpt tpt], [minval-E maxval+E]);
set(h, 'LineStyle', '-', 'Color', [0.6 0.6 0.6]);
set(gca,'ButtonDownFcn',{@click_BSR_ts,handles},'xlim',[0 num_tpts+1]);
setappdata(handles.PrintBSRButton,'bsr_ts',Tseries);


function UpdateTimept(tpt,handles)
% update views
x = getappdata(handles.Sagitalview,'slicenum');
y = getappdata(handles.Coronalview,'slicenum');
z = getappdata(handles.Axialview,'slicenum');
UpdateSlices(x,y,z,handles);

function AxialSlider_Callback(hObject, eventdata, handles)
z = round(get(hObject, 'Value'));
x = round(get(handles.SagitalSlider,'Value'));
y = round(get(handles.CoronalSlider,'Value'));
UpdateSlices(x,y,z,handles);

function SagitalSlider_Callback(hObject, eventdata, handles)
x = round(get(hObject, 'Value'));
z = round(get(handles.AxialSlider,'Value'));
y = round(get(handles.CoronalSlider,'Value'));
UpdateSlices(x,y,z,handles);

function CoronalSlider_Callback(hObject, eventdata, handles)
y = round(get(hObject, 'Value'));
z = round(get(handles.AxialSlider,'Value'));
x = round(get(handles.SagitalSlider,'Value'));
UpdateSlices(x,y,z,handles);

function Timept_Callback(hObject, eventdata, handles)
tpt = round(str2num(get(hObject,'String')));
if tpt<1 || tpt>handles.dims(4)
   prompt={
      ['Invalid Time.'];
      ['Time must be # between: 1 & ',num2str(handles.dims(4))]};
   warndlg(prompt,'Warning:','modal');
end
set(handles.TimeptSlider,'Value',tpt);
UpdateTimept(tpt,handles);

function TimeptSlider_Callback(hObject, eventdata, handles)
tpt = round(get(hObject, 'Value'));
set(handles.Timept,'String',num2str(tpt));
UpdateTimept(tpt,handles);

function click_axial(src,eventdata,handles)%callback for buttonclick
c= ginput(1); 
x = round(c(1,1));
y = round(c(1,2));
z = getappdata(handles.Axialview,'slicenum');
UpdateSlices(x,y,z,handles);

function click_sagital(src,eventdata,handles)
%c = get(gca,'CurrentPoint');
c= ginput(1); 
y = round(c(1,1));
z = round(c(1,2));
x = getappdata(handles.Sagitalview,'slicenum');
UpdateSlices(x,y,z,handles);

function click_coronal(src,eventdata,handles)
%c = get(gca,'CurrentPoint');
c= ginput(1); 
x = round(c(1,1));
z = round(c(1,2));
y = getappdata(handles.Coronalview,'slicenum');
UpdateSlices(x,y,z,handles);

function click_BSR_ts(src,eventdata,handles)
c = ginput(1); 
tpt = round(c(1));
tpt = min(max(1,tpt), handles.dims(4));
set(handles.Timept,'String',num2str(tpt));
set(handles.TimeptSlider,'value',tpt);
UpdateTimept(tpt,handles);

function click_GPAVG_ts(src,eventdata,handles)
c = ginput(1); 
tpt = round(c(1));
tpt = min(max(1,tpt), handles.dims(4));
set(handles.Timept,'String',tpt);
set(handles.TimeptSlider,'value',tpt);
UpdateTimept(tpt,handles);

function Coord_Callback(hObject, eventdata, handles) % user eneters voxel coordinates
coord = round(str2num(get(hObject,'String')));
if (coord(1)<1 || coord(1)>handles.dims(1)) || (coord(2)<1 || coord(2)>handles.dims(2)) || (coord(3)<1 || coord(2)>handles.dims(3))
   prompt={
      ['Invalid Coordinate.'];
      ['X-coord must be # between: 1 & ',num2str(handles.dims(1))];
      ['Y-coord must be # between: 1 & ',num2str(handles.dims(2))];
      ['Z-coord must be # between: 1 & ',num2str(handles.dims(3))]};
   warndlg(prompt,'Warning:','modal');
end
UpdateSlices(coord(1),coord(2),coord(3),handles);

function MNI_coord_Callback(hObject, eventdata, handles)
mni_coord = round(str2num(get(hObject,'String')));
[err, coord] = AFNI_XYZcontinuous2Index(mni_coord, handles.Info, 'LPI', 3);
coord = coord+1; % AFNI indexing convention starts at 0.
UpdateSlices(coord(1),coord(2),coord(3),handles);

function BSR_negthresh_Callback(hObject, eventdata, handles)
thresh = str2num(get(handles.BSR_negthresh,'String'));
if thresh>0, set(handles.BSR_negthresh,'String','0'); end
UpdateBlend(handles);
z = getappdata(handles.Axialview,'slicenum');
x = getappdata(handles.Sagitalview,'slicenum');
y = getappdata(handles.Coronalview,'slicenum');
UpdateSlices(x,y,z,handles);

function BSR_posthresh_Callback(hObject, eventdata, handles)
thresh = str2num(get(handles.BSR_posthresh,'String'));
if thresh<0, set(handles.BSR_posthresh,'String','0'); end
UpdateBlend(handles);
z = getappdata(handles.Axialview,'slicenum');
x = getappdata(handles.Sagitalview,'slicenum');
y = getappdata(handles.Coronalview,'slicenum');
UpdateSlices(x,y,z,handles);

%** Added textbox for clamp range
function TextboxClampRange_Callback(hObject, eventdata, handles)
clamp = str2num(get(handles.TextboxClampRange, 'String'));
if clamp<=0, set(handles.TextboxClampRange, 'String', '10'); end
UpdateBlend(handles);
z = getappdata(handles.Axialview,'slicenum');
x = getappdata(handles.Sagitalview,'slicenum');
y = getappdata(handles.Coronalview,'slicenum');
UpdateSlices(x,y,z,handles);

function SelectMontage_Callback(hObject, eventdata, handles)
contents = cellstr(get(hObject,'String')); % get the whole list as a cell array of strings
montage_type = contents{get(hObject,'Value')}; % get current selection
if strmatch(montage_type,{'Axial','Sagittal','Coronal'})
   TableData = get(handles.MontageViewingTable,'Data'); % get User defined parmaeters for the montage
   figure('NumberTitle','off','name',['Montage ' montage_type]); 
   dims = handles.dims;
   % get blend and subsample according to the montage table parameters
   blend = getappdata(handles.BSR_posthresh,'blend');
   subblend = blend(TableData(1,1):TableData(1,2):TableData(1,3),TableData(2,1):TableData(2,2):TableData(2,3),TableData(3,1):TableData(3,2):TableData(3,3),TableData(4,1):TableData(4,2):TableData(4,3),:); % (x y z t rgb)
  
   % determine axes permuation order 
   switch montage_type
      case 'Axial',    ord = [1 2 3 4]; xlab = 'z slice';
      case 'Sagittal', ord = [2 3 1 4]; xlab = 'x slice';
      case 'Coronal',  ord = [1 3 2 4]; xlab = 'y slice';
   end
 
   % permute subblend in awy that is good for the display
   subblend = permute(subblend,[ord(2) ord(4) ord(1) ord(3) 5]); % 5th dim is rgb
   for d=1:4
      span{d} = TableData(ord(d),1):TableData(ord(d),2):TableData(ord(d),3);
   end
      
   % determine vertical image size (dim1) and horizontal image size (dim2)
   dim1 = numel(span{2}) * numel(span{4}); % rows = ysize*time
   dim2 = numel(span{1}) * numel(span{3}); % columns = xsize*zsize
   
   % reshape subblend and show
   subblend = reshape(subblend,dim1,dim2,3);
   image(subblend,'ButtonDownFcn',{@click_Montage,handles});
   set(gca,'ydir','normal','xtick',([1:numel(span{3})]-0.5)*numel(span{1}),'xticklabel',cellstr(int2str(span{3}'))','ytick',([1:numel(span{4})]-0.5)*numel(span{2}),'yticklabel',cellstr(int2str(span{4}'))');
   ylabel('Time'); 
   xlabel(xlab);
   
   % get current gui (x y z t) coordinate
   z = getappdata(handles.Axialview,'slicenum');
   x = getappdata(handles.Sagitalview,'slicenum');
   y = getappdata(handles.Coronalview,'slicenum');
   t = round(get(handles.TimeptSlider,'value'));
   coord = [x y z t];
   
   % tranaslate this coordinate into this image space
   for d=1:4, subcoord(d) = round((coord(d)-TableData(d,1))/TableData(d,2))+1; end
   subcoord = subcoord(ord);
  
   % plot crosshair and 
   horpos = numel(span{2})*(subcoord(4)-1) + subcoord(2);
   verpos = numel(span{1})*(subcoord(3)-1) + subcoord(1);
   horline = line([1 dim2],[horpos horpos],'color',[0 1 0],'visible','on');
   verline = line([verpos verpos],[1 dim1],'color',[0 1 0],'visible','on');
   
   % record current montage window settings
   setappdata(handles.MontageViewingTable,'span',span);
   setappdata(handles.MontageViewingTable,'ord',ord);
   setappdata(handles.MontageViewingTable,'subcoord',subcoord);
   setappdata(handles.MontageViewingTable,'fig',gcf);
   setappdata(handles.MontageViewingTable,'horline',horline);
   setappdata(handles.MontageViewingTable,'verline',verline);
   setappdata(handles.MontageViewingTable,'horpos',horpos);
   setappdata(handles.MontageViewingTable,'verpos',verpos);
end

function click_Montage(src,eventdata,handles)%callback for buttonclick
% % when user click on the montage figure reverse the caluclation from SelectMontage_Callbac
% % to obtain [x y z t] in the original image space
% c = ginput(1);
% horpos = round(c(1,1));
% verpos = round(c(1,2));
% rspan = getappdata(handles.MontageViewingTable,'span');
% ord = getappdata(handles.MontageViewingTable,'ord');
% TableData = get(handles.MontageViewingTable,'Data'); % get User defined parmaeters for the montage
% revord(ord) = [1:4]; % reverse order; this statement defines new variable revord
% for d=1:4, span{d} = rspan{revord(d)}; end
% subcoord(2) = mod(horpos,numel(span{2})); 
% subcoord(4) = (horpos - subcoord(2))/numel(span{2})+1;
% subcoord(1) = mod(verpos,numel(span{1})); 
% subcoord(3) = (verpos - subcoord(1))/numel(span{1})+1;
% revord(ord) = [1:4]; % reverse order; this statement defines new variable revord
% subcoord 
% ord 
% revord
% %subcoord = subcoord(revord)
% for d=1:4, 
%    coord(d) = round( (subcoord(d)-1)*TableData(d,2) + TableData(d,1) ); 
%    coord(d) = min(max(1,coord(d)),handles.dims(d)); 
% end
% set(handles.Timept,'String',num2str(coord(4)));
% set(handles.TimeptSlider,'value',coord(4));
% 
% % redraw crosshairs
% dim1 = numel(span{2}) * numel(span{4}); % rows = ysize*time
% dim2 = numel(span{1}) * numel(span{3}); % columns = xsize*zsize
% horline = getappdata(handles.MontageViewingTable,'horline');
% verline = getappdata(handles.MontageViewingTable,'verline');
% set(horline,'visible','off');set(verline,'visible','off');
% figure(getappdata(handles.MontageViewingTable,'fig'));
% horline = line([1 dim2],[horpos horpos],'color',[0 1 0],'visible','on');
% verline = line([verpos verpos],[1 dim1],'color',[0 1 0],'visible','on');
% setappdata(handles.MontageViewingTable,'horline',horline);
% setappdata(handles.MontageViewingTable,'verline',verline);
% setappdata(handles.Sagitalview,'slicenum',coord(1));
% setappdata(handles.Coronalview,'slicenum',coord(2));
% setappdata(handles.Axialview,'slicenum',coord(3));
% set(handles.Timept,'String',num2str(coord(4)));
% set(handles.TimeptSlider,'value',coord(4));
% UpdateTimept(coord(4),handles);


function PrintBSRButton_Callback(hObject, eventdata, handles)
val = get(hObject,'Value');
if val==1
   MainVoxTS = getappdata(handles.PrintBSRButton,'bsr_ts')
   set(hObject,'Value',0);
end

%** Also print out timeseries of other datasets.
if ~isempty(handles.OtherImgData)
    x = getappdata(handles.Sagitalview,'slicenum');
    y = getappdata(handles.Coronalview,'slicenum');
    z = getappdata(handles.Axialview,'slicenum');
    Tseries = squeeze(handles.BSR.img(x,y,z,:));

    num_OtherData = length(handles.OtherImgData);
    tmparray = zeros(num_OtherData, length(Tseries));
    for l=1:num_OtherData
        tmparray(l,:) = squeeze(handles.OtherImgData{l}(x,y,z,:));
    end
    tmparray = tmparray';
    
    OtherVoxTS = tmparray
end

function MontageViewingTable_CellEditCallback(hObject, eventdata, handles)
val = round(str2num(eventdata.EditData)); % get cell content as entered by user
inds = eventdata.Indices; % find (i,j) position in the table
% if new value is outside of alowable dims, then change back to previous
% cell content
if ~(val >=1  & val <= handles.dims(inds(1)))
   TableData = get(handles.MontageViewingTable,'Data');
   TableData(inds(1),inds(2)) = eventdata.PreviousData;
   set(handles.MontageViewingTable,'Data',TableData);
end

function MontageViewingTable_CellSelectionCallback(hObject, eventdata, handles)
% hObject    handle to MontageViewingTable (see GCBO)
% eventdata  structure with the following fields (see UITABLE)
%	Indices: row and column indices of the cell(s) currently selecteds
% handles    structure with handles and user data (see GUIDATA)
% NK: note we must have this function even if it doesn't do anything!!!!!





function [cmap,ignore_pts] = set_colormap(max_bsr, min_bsr, pos_bsr_thresh, neg_bsr_thresh)
% written by Jimmy - NK made some changes
%   set the display colormap based on the max/min display values and 
%   the threshold setting.
%
%   The upper colors are coming from the entries of [140:239] of the 
%   255 jet colormap, and the lower colors are from the entries of 
%   [1:100] of the colormap.
%

   range_interval = max_bsr - min_bsr;
   upper_interval = max_bsr - pos_bsr_thresh;
   lower_interval = neg_bsr_thresh - min_bsr; % abs(min_bsr) - abs(neg_bsr_thresh);
   
   %  colormap entries for the upper range values, using the
   %  entries of [140:239] from the 255 jet colormap
   %
   num_upper_colors = 0;
   if (upper_interval > 0)
      num_upper_colors = round(upper_interval / range_interval * 255); 
      cmap_size = round(255 * num_upper_colors/100);
      first_color_idx = round(140 / 255 * cmap_size);
      last_color_idx = first_color_idx + num_upper_colors - 1;
      uppermap = jet(cmap_size);
      upper_colors = uppermap(first_color_idx:last_color_idx,:);
   end;
   
   %  colormap entries for the lower range values, using the 
   %  entries of [1:100] from the 255 jet colormap
   %
   num_lower_colors = 0;
   if (lower_interval > 0)  
      num_lower_colors = round(lower_interval / range_interval * 255); 
      cmap_size = round(255 * num_lower_colors/100);
      first_color_idx = 1;
      last_color_idx = num_lower_colors;
      lowermap = jet(cmap_size);
      lower_colors = lowermap(first_color_idx:last_color_idx,:);
   end;
   
   cmap = zeros(256,3);
   cmap(1:255,:) = jet(255);
  
   ignore_pts = [(num_lower_colors+1):(255-num_upper_colors)];
   
   if (num_lower_colors > 0),
     cmap(1:num_lower_colors,:) = lower_colors;
   end;
   if (num_upper_colors > 0),
     if ~isempty(ignore_pts)
       cmap((ignore_pts(end)+1):255,:) = upper_colors;
     else
       cmap((255-num_upper_colors+1):255,:) = upper_colors;
     end
   end
   if ~isempty(ignore_pts)
     %cmap(ignore_pts,:) = ones(length(ignore_pts),3) * 140/255;
     cmap(ignore_pts,:) = gray(numel(ignore_pts)); 
   end
   cmap(256,:) = [1 1 1];






%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%% Create functions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function BSR_posthresh_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
   set(hObject,'BackgroundColor','white');
end

function BSR_negthresh_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
   set(hObject,'BackgroundColor','white');
end

function MNI_coord_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
   set(hObject,'BackgroundColor','white');
end

function Coord_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
   set(hObject,'BackgroundColor','white');
end

function TimeptSlider_CreateFcn(hObject, eventdata, handles)
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
   set(hObject,'BackgroundColor',[.9 .9 .9]);
end

function Timept_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
   set(hObject,'BackgroundColor','white');
end

function CoronalSlider_CreateFcn(hObject, eventdata, handles)
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
   set(hObject,'BackgroundColor',[.9 .9 .9]);
end

function SagitalSlider_CreateFcn(hObject, eventdata, handles)
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
   set(hObject,'BackgroundColor',[.9 .9 .9]);
end

function AxialSlider_CreateFcn(hObject, eventdata, handles)
if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
   set(hObject,'BackgroundColor',[.9 .9 .9]);
end

function SelectMontage_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end


%** Added clamp range textbox
% --- Executes during object creation, after setting all properties.
function TextboxClampRange_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end